解剖学・骨学実習〜計測・解析編〜
伊藤 毅
12/14/2020
1 はじめに
本文書は、解剖学・骨学実習のうち、計測と解析に関するパートの補足資料です。本実習では、主に霊長類の頭骨標本を材料として、様々な装置を用いた形態計測の方法、3D Slicerを用いたCT画像の処理、Rを用いた幾何学的形態測定について解説します。
間違いがないよう努めてはおりますが、本文書の内容についての責任は負いません。もしお気づきの点等ございましたら、著者までご連絡いただけますと幸いです。
2 計測装置
2.1 装置の種類
- 標本に直接触れて計測する
- ノギス
- MicroScribe（Revware社）：三次元座標データを取得する装置
- ・・・
- 画像データを取得する装置（得られた画像データを用いて様々な計測を行う）
- コンピュータ断層撮影 (CT)
- 磁気共鳴画像（MRI）
- 3Dスキャナ
- デジタルカメラ
- ・・・
ノギス、MicroScribe、一部の3Dスキャナ、デジタルカメラは容易に持ち運びできるので、たとえば博物館を渡り歩いてたくさんのデータを収集するといったことができる。CTやMRIは大掛かりな装置を必要とするが、標本内部の構造を観察することができる。画像データは各種デジタルレポジトリでの整備・公開が進んでおり、世界中の研究者によって有効利用されている。
形態画像データを収録するデジタルレポジトリの例
デジタルカメラで標本を撮像する場合は、視差によって生じる計測誤差に注意する必要がある。これは、標本とレンズの距離を近づけすぎないようにすること、カメラとレンズのセッティングをサンプル間で統一することなどにより、ある程度低減することができるようだ（Mullin & Taylor 2002）。
本実習では、ノギスとMicroScribeの計測を体験する。また、μCT（Bruker社製Skyscan 1275）の撮像を見学する（X線ガラスバッチを着用していない人はCT装置の操作はできない。自分で操作する場合は、事前に所定の講習を受けて登録する必要がある）。
2.2 CTの仕組み
CTは、管球でX線を発生させ、被写体に照射し、反対側に設置された検出器で捉える。このとき、管電圧はX線の強度に、管電流はX線の量に影響する。X線は、被写体によって吸収され、その密度と厚みに応じて減退するので、検出器では被写体のX線透過率の投影分布が取得される。これをあらゆる方向から行い、それらをコンピュータ上で再構成することで、被写体の三次元的構造を計算することができる。以下に、再構成の概念図を示す。 
観測された投影データから未知数Xを推定する。上の単純な例では、線形連立方程式の解を求めることで未知数を特定できる。 \[
X_1+X_2=2\\
X_3+X_4=1\\
X_1+X_3=2\\
X_2+X_4=1\\
X_1+X_4=3
\]
しかし実際には、未知数は巨大になり、また測定値は誤差を含むので、この通りにはいかない。実際のCT装置では、計算負荷の低い、逆投影法を基本とする解析的方法が普及した。以下に、単純逆投影法の概略を示す。 
単純逆投影法は、図に示されるように、被写体が存在しないところにも値が入るので、ボケが生じてしまい実用に適さない。このボケを低減させることを狙いとして、実際のCT装置では、逆投影の前の投影データに（エッジを強調する）フィルタ補正を施す方法（フィルタ補正逆投影法）が一般的に取られてきた。当然、フィルタ（再構成関数）の種類によって再構成画像の印象は大きく変わり、エッジ効果の高いフィルタを適用するほど、空間分解能は高くなるがノイズが増える傾向にある（平尾 2008）。被写体の種類などに応じて、適切な再構成関数を選択する必要がある。
ただし、フィルタ補正を適切に施しても、アーチファクトを完全に消すことはできない。近年では、コンピュータの性能の向上に伴って、代数的・統計的再構成方法が採用されるようになり、ノイズやアーチファクトの低減といった、従来の解析的方法の欠点の解決が期待されている。
2.3 参考・引用文献
- 田中敏幸. (2019) X線CTの撮影原理. 計測と制御 58: 509-513.
- 齊藤泰男. (2009) X線CT 第1回：画像形成の原理，装置開発の現状. Medical Imaging Technology 27: 200-204.
- 平尾芳樹. (2008) 医療用X線CT技術の系統化調査報告. 国立科学博物館技術の系統化調査報告 12: 85-160.
- 辻岡勝美, 井田義宏. (2002) CT 検査の実際 (1): スキャンパラメータの設定 (< 基礎連続講座> CT 講座). 日本放射線技術学会雑誌 58: 1456-1460.
- Mullin SK, Taylor PJ (2002) The effects of parallax on geometric morphometric data. Comput Biol Med 32: 455–464
3 CT画像解析
3.1 CT画像について
再構成されたデータは、通常、TIFFなどの一般的な画像形式もしくはDICOMという医用画像の標準規格をみたす形式で、2D画像のスタックという形で出力される。各ピクセルはX線吸収係数を換算したCT値を持ち、通常、水が0、空気が-1000、となるようにキャリブレーションされている（Hounsfield unit: HU）。CT値が大きいほど、密度が高く、X線吸収係数が大きいことを表す。
3.2 CT画像の表示
以下に示すように、さまざまなソフトウェアがある。
- Amira - Avizo（有償）
- Horos（フリー、macのみ）
- 3D Slicer（フリー、クロスプラットフォーム）
- Fiji - ImageJ（フリー、クロスプラットフォーム）
- ・・・
本実習では3D Slicer 4.11を利用する。
Digital Morphology Museum, KUPRIからダウンロードした以下のデータを使用する。
| PRICT | ID | Species | Abbreviation |
|---|---|---|---|
| 504 | WPK-Mimi | Pan troglodytes | Pt |
| 949 | PRI-673 | Cebus albifrons | Ca |
| 1055 | PRI-4536 | Presbytis melalophos | Pm |
| 1095 | PRI-6093 | Varecia variegata | Vv |
| 1267 | PRI-6494 | Macaca fuscata | Mf |
| 1287 | PRI-2516 | Hylobates lar | Hl |
dicomの入ったフォルダの選択 
- 2D Slicerの“Load DICOM Data”を開く
dicomの読み込み 
- “Import DICOM Data”で、dicomデータの入ったフォルダを選択する。
- データを選択する。
- “Load”で読みこむ。
2D表示 
- スクロールバーを動かすことで、スライス間を移動できる
- ここを選択してから画面上をドラッグすると、Window Level（明るさ）とWindow Width（コントラスト）を調整できる。縦方向のドラッグでLevel、横方向のドラッグでWidthを調整する。
3D表示（ボリュームレンダリング） 
- “Volume Rendering”ツールを選択する
- 眼のアイコンをクリックして開く
- “Preset”から画像や被写体の種類に応じてカラーを選択する
- ここをクリックすると3D画像が中心位置に戻される（これをしないとなぜか3Dが表示されない…）
- “Shift”のスクロールバーで閾値の変更ができる
3.3 3Dサーフェスモデルの作成
3Dサーフェスモデルを作成するためには、まずセグメンテーション（画像の2値化）を行う必要がある。
本実習では、もっとも簡単なグローバル閾値法によるセグメンテーションを行う。今回、閾値は、見た目で判断して適当に設定するが、厳密な精度が要求される場合はhalf-maximum height (HMH)法（Spoor et al. 1993; Coleman and Colbert 2007）により求めるのが望ましい。たとえば、Coleman and Colbert (2007)では、ランダムに選択した10箇所においてHMHを算出し、それらの平均値をグローバル閾値としている。
HMHは、境界領域におけるCT値の最大値と最小値の平均値である。以下に例を示す。
空気と骨の境界領域に引いた線分上で、CT値の分布を示す。 
ここで、最大値は1755、最小値は-1025なので、HMHはそれらの平均値として365となる。
HMH（365）をグローバル閾値をとしてセグメンテーションすると、以下の図のようになる。 
仮に、閾値を-900とすると、以下の図のようになる。HMHのそれよりも骨領域が膨らんでいることがわかる。 
一方、閾値を900とすると、以下の図のように、骨領域が細くなる。 
このように、閾値は、高すぎると実際よりも痩せてしまったり穴が空いてしまったりし、逆に低すぎると実際よりも膨らんでしまったりアーチファクトが映り込んでしまったりする。このような性質を持つため、CT値のバラツキが大きい領域をセグメンテーションしようとする場合、グローバル閾値法は適さない。グローバル閾値法のほか、以下に示すような様々なセグメンテーションの方法が開発されており、データの性質や求められる精度に応じて使い分けるとよい。各手法の適用例や精度評価については、Engelkes (2020)、Rathnayaka et al. (2011)、Scherf & Tilgner (2009)、Ito (2019)などに解説がある。
- 手動
- ローカル閾値
- エッジ検出・分水嶺
- レジストレーション
- 深層学習
- ランダムウォーク
- ・・・
グローバル閾値法によるセグメンテーション 
- “Editor”モジュール（“Legacy”の下層にある）を選択する。color tableを選択するよう指示が出るが、デフォルトの設定のまま“OK”をクリックする。
- “Master Volume”が目的のデータになっていることを確認する
- “ThresholdEffect”のアイコンをクリックする
- “Label”の色を選択する。今回は骨なので“bone”にする
- “Threshold Range”を設定する（グローバル閾値の設定）。
- “Apply”をクリックする
サーフェスレンダリング 
- “MakeModelEffect”のアイコンをクリックする
- “Apply”をクリックする
- ここをクリックすると3D画像が中心位置に戻される（これをしないとなぜか3Dが表示されない…）
サーフェスモデルの保存 
- “Save”をクリックし、データの選択、必要に応じて名前の変更、パスの変更を行い保存する。保存形式は、stlかplyでいいだろう。
3.4 ランドマークデータの取得
ランドマークデータを取得するソフトウェアには、例えば以下のものがある。セミランドマークについては後述する。
- 3D Slicer (3D、セミランドマーク可)
- TpsDig2 (Windowsのみ、2Dのみ、セミランドマーク可)
- geomorph (Rパッケージ、2D3D、セミランドマーク可)
- Checkpoint (有償、Windowsのみ、3D、セミランドマーク可)
- Amira - Avizo (有償、2D3D)
本実習では再び3D Slicer 4.11を使用し、以下のランドマークのデータを取得する。
| Order | Landmark | Type |
|---|---|---|
| 1 | Prosthion | 2 |
| 2 | Glabella | 3 |
| 3 | Inion | 2 |
| 4 | Basion | 2 |
| 5 | Zygomaxillare (left) | 2 |
| 6 | Zygomaxillare (right) | 2 |
Bookstein (1991)は、ランドマークを3つのタイプに分類している。
- Type 1: discrete juxtapositions of tissues
- Type 2: maxima of curvature or other local morphogenetic processes
- Type 3: extremal points
幾何学的形態測定（次節で解説）では、ランドマークが相同であることを前提とする。したがって、確実に相同性を担保できるタイプ1のランドマークがもっとも望ましく、そうではないタイプ3のランドマークはあまり望ましくない。タイプ2のランドマークはそれらの中間的性質をもつ。
ランドマーキング 
- 表示モードを“3D only”にすると見やすい
- “Markups”モジュールを選択する
- “Fiducial”を選択。マウスのカーソルを目的のランドマークの位置に移動させクリックするとその位置の座標データを取得できる。
- “Display”のスクロールバーでランドマークのサイズや文字サイズを変更できる
ランドマークデータの保存 
- 3Dサーフェスモデルを保存した時と同様に、“Save”をクリックし、データの選択、必要に応じて名前の変更、パスの変更を行い保存する。保存形式は、座標データのみであれば、fcsvでよいだろう（警告メッセージが出るがOKとする）。本実習では、fcsvで保存する。表示形式も保存しておく必要がある場合は、json形式で保存する。
fscvファイルの中を見てみよう。
## # Markups fiducial file version = 4.11
## # CoordinateSystem = LPS
## # columns = id,x,y,z,ow,ox,oy,oz,vis,sel,lock,label,desc,associatedNodeID
## 1,-2.0602643489837646,-87.3006591796875,-1103.945068359375,0,0,0,1,1,1,0,F-1,,
## 2,-0.43642905354499817,-30.19694709777832,-1016.7828979492188,0,0,0,1,1,1,0,F-2,,
## 3,2.1217479705810547,100.34082794189453,-1020.6141967773438,0,0,0,1,1,1,0,F-3,,
## 4,3.224494457244873,62.182159423828125,-1080.2637939453125,0,0,0,1,1,1,0,F-4,,
## 5,46.13127517700195,-23.04610824584961,-1078.1734619140625,0,0,0,1,1,1,0,F-5,,
## 6,-42.983951568603516,-23.148298263549805,-1081.05810546875,0,0,0,1,1,1,0,F-6,,はじめの3行はコメントで、4行目から取得した順番に座標データが並ぶ。1列目は番号で、2〜3列目にXYZ座標データがミリ単位で記録されている。
次節では、このデータをRに読み込んで、幾何学的形態測定による解析を行う。
3.5 参考・引用文献
- 塚越伸介. (2009) X線CT第2回:CT画像の基本と画像表示. Medical Imaging Technology 27: 258-262.
- Spoor C.F., Zonneveld F.W., & Macho G.A. (1993) Linear measurements of cortical bone and dental enamel by computed tomography: Applications and problems. American Journal of Physical Anthropology 91: 469–484.
- Coleman M. and Colbert M. (2009) Technical note: CT thresholding protocols for taking measurements on three‐dimensional models. American Journal of Physical Anthropology 133: 723-725.
- Engelkes K. (2020) Accuracy of bone segmentation and surface generation strategies analyzed by using synthetic CT volumes. https://doi.org/10.13140/RG.2.2.18389.86247
- Rathnayaka K., Sahama T., Schuetz M.A., & Schmutz B. (2011) Effects of CT image segmentation methods on the accuracy of long bone 3D reconstructions. Medical Engineering & Physics 33: 226–233.
- Scherf H. & Tilgner R. (2009) A new high-resolution computed tomog- raphy (CT) segmentation method for trabecular bone architectural analysis. American Journal of Physical Anthropology 140: 39–51.
- Ito T. (2019) Effects of different segmentation methods on geometric morphometric data collection from primate skulls. Methods in Ecology and Evolution 10: 1972-1984.
- https://www.youtube.com/watch?v=MKLWzD0PiIc
4 幾何学的形態測定
4.1 幾何学的形態測定とは
幾何学的形態測定（geometric morphometrics）では、古典的な形態測定法のように距離や角度に落とし込むことなく、座標データそのものを扱う。幾何学的形態測定では、拡大・縮小（scaling）、 回転（rotation）、移動（translation/centering）といった座標変換を施すことで、サンプル間の大きさ、方向、位置の違いを調整する。この処理は、Procrustes整列（Procrustes superimposition）、Procrustes fitting、Generalized Procrustes Analysis (GPA)などと呼ばれる。このとき、大きさ（サイズ）は重心サイズ（Centroid size；重心から各標識点までの二乗距離の総和の平方根）で表される。Procrustes整列後のデータを形状（シェイプ）データとして利用する。
シェイプはサイズの成分を含まない（重心サイズが1に揃えられている）。ただしこれは、シェイプがサイズと無相関であるということではない。サイズと無相関のシェイプ成分を扱いたい場合は、シェイプをサイズに回帰させ、その残差をとる。サイズと相関するシェイプはアロメトリー成分という。
幾何学的形態測定については、Bookstein (1991)、Zelditch et al. (2012)、Mitteroecker & Gunz (2009)、Slice (2007)、Adams et al. (2004)などの参考文献がある。
4.2 基本の解析
今回の実習では、簡単のため6サンプルのみを解析に用いるが、ランドマーク数に対してサンプルサイズが少なすぎると統計的な問題となる場合があるので、本来は望ましくない。実際の研究では、十分なサンプルサイズを確保するようにしたい。
手始めに、Pan troglodytesのランドマークデータを読み込んでみる。
library(tidyverse)
read_csv("../landmark/Pan_troglodytes.fcsv",
skip = 3, # 最初の3行のコメント行をスキップ
col_names = FALSE)[,2:4] # 座標データが含まれる2,3,4列を抽出## # A tibble: 6 x 3
## X2 X3 X4
## <dbl> <dbl> <dbl>
## 1 -2.06 -87.3 -1104.
## 2 -0.436 -30.2 -1017.
## 3 2.12 100. -1021.
## 4 3.22 62.2 -1080.
## 5 46.1 -23.0 -1078.
## 6 -43.0 -23.1 -1081.次に、関数を作成し、指定したフォルダに含まれるすべてのデータを一挙に読み込み配列（array）に格納する。
read_fcsv <- function(path, k = 3, p = 6){
list_data <- lapply(
list.files(path, full.names = TRUE),
function(x){read_csv(x, skip = 3, col_names = FALSE)}[,2:4]%>%as.matrix)
# 配列に格納する
array_data <- array(as.numeric(unlist(list_data)), dim=c(p, k, length(list_data)))
# 名前をつける
names <- tools::file_path_sans_ext(list.files(path, pattern = "*.csv"))
dimnames(array_data)[[3]] <- names
return (array_data)
}
raw_data <- read_fcsv("../landmark")
dimnames(raw_data)[[3]] <- c("Ca", "Hl", "Mf", "Pt", "Pm", "Vv") # 種名が長いので略称に
raw_data## , , Ca
##
## [,1] [,2] [,3]
## [1,] -0.1149868 -35.12373 -1054.850
## [2,] -0.6637123 -17.67938 -1023.330
## [3,] -0.1195280 53.23666 -1045.329
## [4,] -0.1068767 23.19316 -1060.824
## [5,] 19.9068317 -13.63381 -1053.072
## [6,] -20.5708523 -13.87484 -1053.465
##
## , , Hl
##
## [,1] [,2] [,3]
## [1,] 1.4124396 49.33804 -1338.383
## [2,] -0.2931501 30.91640 -1305.041
## [3,] -0.8856894 -54.35939 -1319.973
## [4,] 2.1998587 -23.15752 -1341.715
## [5,] -23.5790615 21.08697 -1340.217
## [6,] 23.8187637 21.20558 -1338.248
##
## , , Mf
##
## [,1] [,2] [,3]
## [1,] -0.8745944 -51.14996 -1061.536
## [2,] -0.9600250 -17.43713 -1020.646
## [3,] -1.9526790 65.18864 -1053.153
## [4,] -1.5969342 28.84805 -1069.347
## [5,] 27.5209427 -15.30799 -1064.343
## [6,] -27.5567760 -15.59646 -1066.165
##
## , , Pt
##
## [,1] [,2] [,3]
## [1,] -2.0602643 -87.30066 -1103.945
## [2,] -0.4364291 -30.19695 -1016.783
## [3,] 2.1217480 100.34083 -1020.614
## [4,] 3.2244945 62.18216 -1080.264
## [5,] 46.1312752 -23.04611 -1078.173
## [6,] -42.9839516 -23.14830 -1081.058
##
## , , Pm
##
## [,1] [,2] [,3]
## [1,] -2.2882662 -38.44722 -1014.9694
## [2,] -0.9373286 -16.93668 -985.1452
## [3,] -0.3596495 54.95614 -1010.1675
## [4,] -0.9612941 25.87814 -1026.7108
## [5,] 24.1622162 -19.22473 -1017.1849
## [6,] -27.1201286 -17.04971 -1016.9243
##
## , , Vv
##
## [,1] [,2] [,3]
## [1,] 0.1620082 -52.117371 -1306.246
## [2,] -0.9636781 -15.641431 -1282.679
## [3,] -2.3839674 53.863892 -1289.016
## [4,] -0.2980522 41.737179 -1307.885
## [5,] 20.6963501 -3.194123 -1302.389
## [6,] -21.0222626 -3.935471 -1304.899読み込んだ座標データを表示してみる。
library(geomorph)
plotAllSpecimens(raw_data)
CT撮像時に標本の位置を揃えていないため、とくにZ方向の位置がバラバラである。方向は概ね揃えてあるが、厳密には揃っていない。当然、大きさもバラバラである。小さな黒い点は平均形状を示す。
そこで、Procrustes整列を行う。
gpa_data <- gpagen(raw_data)##
|
| | 0%
|
|============== | 20%
|
|============================ | 40%
|
|======================================================================| 100%plotAllSpecimens(gpa_data$coords)
サンプル間の大きさ、方向、位置の違いが調整され、整列されていることがわかる。
整列後のデータに対して、主成分分析を実行する。第一主成分と第三主成分の散布図を描く。
pca_result <- gm.prcomp(gpa_data$coords)
library(ggplot2)
pca_result$x %>%
as_tibble() %>%
mutate(Species = dimnames(raw_data)[[3]]) %>%
ggplot(aes(x = Comp1, y = Comp2, label = Species)) +
geom_point() +
geom_text(hjust = 0, nudge_x = 0.005)
ちなみに、Rのデフォルト関数（prcomp）を使用して、prcomp(two.d.array(gpa_data$coords))としても同一の結果が得られる。
各主成分の寄与率（proportion of variance）を表示する。
tibble(
PC = 1:6,
Proportion_of_variance = pca_result$sdev^2/sum(pca_result$sdev^2),
Cummurative_proportion = cumsum(pca_result$sdev^2/sum(pca_result$sdev^2))
) %>%
mutate(across(2:3, round, 2))## # A tibble: 6 x 3
## PC Proportion_of_variance Cummurative_proportion
## <int> <dbl> <dbl>
## 1 1 0.7 0.7
## 2 2 0.18 0.88
## 3 3 0.07 0.95
## 4 4 0.04 0.99
## 5 5 0.01 1
## 6 6 0 1第一主成分と第二主成分で、全主成分のうち88%の分散を説明している。
第一主成分の形状変化を可視化する。平均形状から第一主成分の最大値への変化を示す。
plotRefToTarget(gpa_data$consensus, pca_result$shapes$shapes.comp1$max, method = "vector", label = TRUE)
第一主成分が大きいほど、全体的に細長い傾向にあるようだ。
続いて、第二主成分。
plotRefToTarget(gpa_data$consensus, pca_result$shapes$shapes.comp2$max, method = "vector", label = TRUE)
第二主成分が大きいほど、相対的に顔面高が大きいようだ。
主成分得点は、回帰分析や分散分析など、さまざまな一般的な統計分析に利用できる。
4.3 系統解析
形状データの距離行列を用いて、系統関係を推定してみよう。
手始めに、CaとHlの間のユークリッド距離を計算してみる。 \[d = \sqrt{\sum(X_i - Y_i)^2} \]
sqrt(
sum(
(pca_result$x[1, ] -
pca_result$x[2, ]
)^2
)
)## [1] 0.0613446これは、dist関数を用いることで、一挙に計算できる。
dist(pca_result$x)## Ca Hl Mf Pt Pm
## Hl 0.06134460
## Mf 0.10520024 0.09636327
## Pt 0.17410217 0.17332358 0.13863467
## Pm 0.08351498 0.09022352 0.12018503 0.17620024
## Vv 0.28935364 0.25498915 0.25883128 0.22524408 0.29382293これ（すべての主成分の得点のユークリッド距離）は、Procrustes距離と一致する。
gpa_data$procD## Ca Hl Mf Pt Pm
## Hl 0.06134460
## Mf 0.10520024 0.09636327
## Pt 0.17410217 0.17332358 0.13863467
## Pm 0.08351498 0.09022352 0.12018503 0.17620024
## Vv 0.28935364 0.25498915 0.25883128 0.22524408 0.29382293これで距離行列の準備ができた。では、まずは、もっとも単純な方法の一つである、UPGMA（Unweighted pair group method with arithmetic mean）法で系統を推定してみよう。
このうち、最小の距離をもつ組み合わせは、CaとHl（0.0613446）である。これらを一つにまとめ、クラスタとし、距離行列をアップデートする。ここで、各要素とクラスタとの距離は、そのクラスタに含まれる全ての要素との距離の平均とする。
## Mf Pt Pm Vv
## Pt 0.13863467
## Pm 0.12018503 0.17620024
## Vv 0.25883128 0.22524408 0.29382293
## (Ca,Hl) 0.10078175 0.17371287 0.08686925 0.27217139枝長は距離の半分として0.0306723となる。 
アップデートされた距離行列の中で、最小の距離をもつ組み合わせは、Pmと(Ca,Hl)（0.10078175）である。これらをまとめると次のようになる。
## Mf Pt Vv
## Pt 0.1386347
## Vv 0.2588313 0.2252441
## (Pm,(Ca,Hl)) 0.1072495 0.1745420 0.2793886
以下同様に、最小の距離をもつものから、クラスタにまとめ上げていく。
## Pt Vv
## Vv 0.2252441
## (Mf,(Pm,(Ca,Hl))) 0.1655652 0.2742493
UPGMAは、hclust関数を使って計算することもできる。
library(ape)
plot(
as.phylo(
hclust(dist(pca_result$x), method = "average")
)
)
axisPhylo()
UPGMAは進化速度の一定性を仮定しており、有根系統樹が得られる。
一方、進化速度一定の仮定を要しない方法としては、近隣結合法（neighbor joining）がある。
nj_morph <- nj(dist(pca_result$x))
plot(nj_morph ,type="unrooted")
Vvを外群とした時の有根系統樹（頭蓋形状）
nj_morph_rooted <- root(nj_morph, "Vv")
plot(nj_morph_rooted)
axisPhylo()
UPGMAと概ね同じ分岐関係が得られるが、Hl、Ca、Pmの分岐関係は異なっている。
では、頭蓋形状のデータから推定した系統は、系統関係を正しく推定できているだろうか。 DNAの塩基配列に基づく系統推定の結果と比較してみよう。 10kTreesからアライメント済のデータ [Final dataset (NEXUS format)] をダウンロードして、Rに読み込み、近隣結合法により系統を推定する。
# 10kTreesのnexus dataの読み込み
nexus_10ktrees <- read.nexus.data("../10kTrees/10kTrees_finalFile_version3.nex")
# 形態データと同じ6種のデータを抽出
nexus_sub <- nexus_10ktrees[c("Cebus_albifrons", "Macaca_fuscata", "Pan_troglodytes_verus", "Presbytis_melalophos", "Varecia_variegata_variegata", "Hylobates_lar")]
# DNAbin形式に変換
dnabin_sub <- nexus2DNAbin(nexus_sub)
# 名前を略称に
names(dnabin_sub) <- c("Ca", "Mf", "Pt", "Pm", "Vv", "Hl")
nj_dna <- nj(dist.dna(dnabin_sub))
plot(nj_dna,type="unrooted")
Vvを外群とした時の有根系統樹（DNA塩基配列）
nj_dna_rooted <- root(nj_dna, "Vv")
plot(nj_dna_rooted)
axisPhylo()
DNA塩基配列に基づく系統樹では、旧世界ザル類のMfとPm、ヒト上科のPtとHlがそれぞれクラスターを形成し、その外側に新世界ザル類のCaが来ており、枝長はともかく分岐関係は正しく推定できているように見える。一方、形状に基づく系統樹は、まったく間違った分岐関係になっている。（少なくとも今回のデータで表される）頭蓋形状の変異は、系統関係を反映しないようだ。
幾何学的形態測定のデータに基づく系統推定の方法論については、未だ活発に議論が続いている。よくある落とし穴は、一部の主成分の得点を系統解析の形質として用いるというもので、問題があるようだ（Adams et al. 2011）。本実習で紹介した、Procrustes距離行列に基づくUPGMAや近隣結合法のほか、最節約法による推定方法も開発されている（Catalano et al. 2011）。
最後に、形状の主成分得点の散布図に、DNA塩基配列に基づく系統樹をプロットしてみよう。
pca_result <- gm.prcomp(gpa_data$coords, phy = nj_dna)
plot(pca_result, phylo = TRUE)
これはphylomorphospaceと呼ばれる解析手法で、主成分で表される形状空間において、系統発生の過程でどのような変化があったのかが可視化される。これは単に主成分空間に後から系統樹をプロットしたもので、主成分自体は先に示したものと同一である。
4.4 計測誤差
（準備中）
4.5 欠損値の補完
（準備中）
4.6 セミランドマーク
（準備中）
- Procrustes距離最小化法
- 屈曲エネルギー最小化法
4.7 参考・引用文献
- Bookstein F.L. (1991) Morphometric Tools for Landmark Data: Geometry and Biology. Cambridge University Press, Cambridge.
- Zelditch M.L., Swiderski D.L., Sheets H.D. (2012) Geometric Morphometrics for Biologists: A Primer. Academic Press, New York.
- Mitteroecker P. & Gunz P. (2009) Advances in geometric morphometrics. Evolutionary Biology 36: 235-247.
- Slice D.E. (2007) Geometric morphometrics. Annual Review of Anthropology 36: 261-281.
- Adams D.C., Rohlf F.J., Slice D.E. (2004) Geometric morphometrics: ten years of progress following the ‘revolution’. Hystrix: Italian Journal of Zoology 71: 5-16.
- Arnold C., Matthews L.J., Nunn C.L. (2010) The 10kTrees Website : A New Online Resource for Primate Phylogeny. Evolutionary Anthropology 19: 114-118.
- Adams D.C., Cardini A., Monteiro L.R., O’Higgins P., Rohlf F.J. (2011) Morphometrics and phylogenetics: Principal components of shape from cranial modules are neither appropriate nor effective cladistic characters. Journal of Human Evolution 60: 240-243.
- Catalano S.A., Goloboff P.A., Giannini N.P. (2010) Phylogenetic morphometrics (I): The use of landmark data in a phylogenetic framework. Cladistics 26: 539-549.